Phase diagram and structure of colloid-polymer mixtures confined between walls 
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The influence of confinement, due to flat parallel structureless walls, on phase separation in 
colloid-polymer mixtures, is investigated by means of grand-canonical Monte Carlo simulations. 
Ultra-thin films, with thicknesses between D = 3 — 10 colloid diameters, are studied. The Asakura- 
Oosawa model [J. Chem. Phys. 22, 1255 (1954)] is used to describe the particle interactions. To 
simulate efficiently, a "cluster move" [J. Chem. Phys. 121, 3253 (2004)] is used in conjunction 
with successive umbrella sampling [J. Chem. Phys. 120, 10925 (2004)]. These techniques, when 
combined with finite size scaling, enable an accurate determination of the unmixing binodal. Our 
results show that the critical behavior of the confined mixture is described by "effective" critical 
exponents, which gradually develop from values near those of the two-dimensional Ising model, to 
those of the three-dimensional Ising model, as D increases. The scaling predictions of Fisher and 
Nakanishi [J. Chem. Phys. 75, 5875 (1981)] for the shift of the critical point are compatible with our 
simulation results. Surprisingly, however, the colloid packing fraction at criticality approaches its 
bulk {D — > oo) value non-monotonically, as D is increased. Far from the critical point, our results 
are compatible with the simple Kelvin equation, implying a shift of order 1/D in the coexistence 
colloid chemical potential. We also present density profiles and pair distribution functions for a 
number of state points on the binodal, and the influence of the colloid-wall interaction is studied. 

PACS numbers: 64.70Fx, 64.60Fr, 05.70.Jk 



I. INTRODUCTION 



Confinement of fluids in nanoscopic capillaries is a 
problem that has received longstanding attention. An in- 
teresting interplay occurs between surface effects at the 
confining walls, such as wetting or drying P,|3,0iBj and 
finite size effects due to the finite capillary width. This 
interplay leads to a host of intriguing phenomena, such as 
capillary condensation or evaporation |^, V, "P] , and a 
crossover in critical behavior from three-dimensional (3d ) 
to two-dimensional (2d) Ising character 0, |3, ITol Hi! ] . 
Apart from fundamental theoretical importance, under- 
standing the structure and phase behavior of confined 
fluids is also a prerequisite for targeted applications in 
microfluidic and nanofluidic devices, which are becoming 
increasingly more relevant LaLJ,UJ,aJ. 

However, for fluids consisting of small molecules, com- 
plications such as the atomistic corrugation of the confin- 
ing walls [iSl ]. impurity atoms at the walls, roughness of 
the walls due to surface steps, dislocations in the crystal 
structure, and so forth, may have a profound effect on 
the phenomena mentioned above. This is especially true 
when the width of the slit pore is of the order of a few 
nanometers, in which case a quantitative understanding 
becomes difficult. In this respect, colloidal dispersions, 
containing colloidal particles with diameters in the mi- 
crometer range, possess certain advantages. It then be- 
comes possible to prepare slit pores which are only a few 
particle diameters wide, and yet have walls that are es- 
sentially fiat on the size scale of the particles. In addi- 
tion, interactions between colloidal particles can be well 
tuned [T^lrRITsj . Particularly promising systems of this 



kind are colloid-polymer mixtures, since both the bulk 
phase behavior, and the interfaces separating the colloid- 
rich and polymer-rich phases, can be studied experimen- 
tally in detafi [l9 l [2 ^j 2i y 2al . Moreover, the Asakura- 
Oosawa (AO) modeljillil provides a simple theoretical 
framework, which seems to capture the essential features 
of such phase-separating systems, and is moreover well 
suited for computer simulations IM El IHIll El 113 ■ 

The aim of this work is to provide precise predictions 
for the phase diagram and the structure of the AO model 
confined between structureless hard walls. We consider 
film thicknesses ranging from D = 3 to I? = 10 colloid 
diameters. Our work complements earlier work E^jEIIj 
based on Gibbs ensemble Monte Carlo, which focused 
on the non-critical regime of the phase diagram. In this 
work, also the critical regime is considered, in order to 
compare to theoretical predictions for the shift in critical 
parameters as function of the film thickness H . We 
consider a colloid-to-polymer size ratio q = ap/ac = 0.8, 
with (Tc the colloid diameter and dp the polymer diameter 
of gyration, since for this particular size ratio, accurate 
information on the bulk and intcrfacial properties is avail- 
able [2I m . We also explore the effect of an additional 
weak repulsion between the walls and the colloids (in ad- 
dition to the hard interaction already present). Under- 
standing the combined effect of film thickness and wall- 
colloid interaction, is crucial in providing guidance for the 
interpretation of possible experiments. At this point, we 
are aware of only one experiment on capillary condensa- 
tion of a colloid-polymer mixture (in a wedge formed by 
glass plates ) , but we hope that the present study will 
encourage further experiments. In related previous work 
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carried out by us, in which a film thickness D = 5 and 
purely hard walls were considered, the crossover from 3d 
to 2d Ising critical behavior was already discussed 0- 

The outline of this paper is as follows. In Sec. II, 
we summarize the most important theoretical predictions 
relevant for the interpretation of our results. In Sec. Ill, 
the AO model is introduced, and details on the simula- 
tion technique are provided. In Sec. IV, we investigate 
the phase behavior and the structural properties of the 
AO model as the film thickness D is varied, while Sec. V 
considers the influence of the colloid-wall interaction. Fi- 
nally, in Sec. VI, we summarize our main conclusions. 

II. THEORETICAL BACKGROUND 

We consider a system with an one-component order pa- 
rameter near its critical point. Examples of such systems 
are fluids near their vapor-liquid critical point, magnetic 
systems (such as the simple Ising model), and binary 
mixtures (such as colloid-polymer mixtures) near their 
critical point of unmixing. As is well known, all these 
systems in 3d bulk belong to the 3d Ising universality 
class In particular, critical exponents such as the 

exponent v (which characterizes the growth of the order 
parameter correlation length ^ near the critical point), 
or the order parameter exponent P (which characterizes 
the shape of the binodal between the unmixed phases 
in a binary mixture) have nontrivial values v « 0.630, 
/3 « 0.326 1^1 El 0, rather than the classical mean- 
field values fMF = 1/2, Pmf = 1/2. 

When these systems are confined between two identi- 
cal flat planar structureless walls a distance D apart, two 
distinct physical phenomena affecting the phase transi- 
tion are expected to occur: 

(i) The growth of the critical correlations in the direc- 
tion perpendicular to the walls is limited by the fi- 
nite thickness of the film, while critical correlations 
can still grow further and unlimited in the two di- 
rections parallel to the film. Therefore, a crossover 
from 3d critical behavior toward 2d critical behav- 
ior (for which the exponents mentioned above take 
the values i^ -y = 1 and P2 ~ 1/8, assuming 2d Ising 
universality lU ) is expected. This crossover is 
predicted to occur when the temperaturc-like vari- 
able has approached a relative distance t from the 
critical point in the bulk such that ^ cx t^'^ is of the 
same order as D, implying a crossover distance of 
order icross D^^/'^ . This dimensional crossover 
in critical behavior also implies an additional en- 
hancement of fluctuations (in lower dimensionality 
stronger fluctuations occur). This leads to a de- 
pression of the critical point. Consequently, one 
expects a shift of the critical point of the same or- 
der, Uift K 

(ii) If there is no special symmetry between the coex- 
isting phases, one must expect that there is an (en- 



ergetic and/or entropic) preference of the walls for 
one of the coexisting phases in comparison to the 
other. This phenomenon is also well known to oc- 
cur away from the critical point. For example, hy- 
drophilic walls of a thin slit capillary are known 
to lead to "capillary condensation" of an under- 
saturated vapor in the capillary |36l |. Of course, 
for hydrophobic rather than hydrophilic walls, also 
the opposite effect ( "capillary evaporation" ) may 
occur. 

Clearly, a quantitative prediction of the magnitude, 
and sometimes even the sign, of these effects, requires de- 
tailed knowledge of the forces between the walls and the 
particles confined by them. However, since phase equi- 
librium is always shifted due to surface corrections to the 
relevant thermodynamic potentials, we expect a shift in 
the coexistence chemical potential in the thin film, rela- 
tive to its value in the bulk, of order 1/D. This is the 
so-called "Kelvin equation" |j,[a,Sl3i which is expected 
to be valid away from the critical point. In contrast, the 
corresponding shift in the chemical potential at critical- 
ity is more subtle. This becomes most transparent for 
the case of an Ising magnet (or, equivalently, the lattice 
gas model), where in the bulk there is a symmetry be- 
tween the coexisting "spin up" and "spin down" phases. 
However, in the thin film, this symm etry may be broken 
by a surface magnetic field i?i H i, IHl HHll . While 
phase coexistence in the bulk occurs for bulk magnetic 
field H = 0, coexistence in the thin film now requires a 
nonzero bulk field (with a sign opposite to that of the 
surface field Hi). 

It has been shown that in the presence of the surface 
field, the singular part of the free en ergy per spin in this 
Ising model should scale as |1 15II11II3I l38|| 

UnsiD,T,H,Hi)^ (1) 

with a « 0.110 [33. [3^ I33 the critical exponent of the 
specific heat, A w 1.56 |32l IssL l33 | the so-called "gap 
exponent" which characterizes the bulk equation of state, 
and Ai « 0.47 ^38, 39, 40, 4^ its surface analogue. The 
scaling functions /-t, where the signs refer to the sign 
of i, are discussed in more detail in Refs. Here, 
we only infer from Eq. ^ that the shift in the critical 
temperature should read as 

Uift = AT,/T,{^) (2) 
= [rc(i^,i/i)-rc(oo)]/Te(oo) 

and similarly for the shift in the bulk magnetic field re- 
quired to establish coexistence again 

AH, = H,{D,Hi) (3) 
= -D-^/^Y^iHiD^'/") 
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Here, Xc and are again scaling functions, and rc(oo) 
denotes the critical temperature in the bulk [D — > oo) 
system. Note that for small Hi , we expect that Xc tends 
to a constant, while Yc becomes a linear function of its 
argument. This limit is relevant in the present case, since 
we deal with rather small film thicknesses D. 

Considering now mixtures of colloids (c) and poly- 
mers (p) in the grand-canonical ensemble, chemical po- 
tentials, and /ip, of colloids and polymers, respec- 
tively, as well as temperature T and system volume V 
become the relevant independent thermodynamic vari- 
ables. In the framework of the AO model [23| . 
the dependence on temperature only enters via the fu- 
gacities = exp(/Xc/fcB'?^)j and Zp = exp(/Xp/fcBT). 
There is no other explicit temperature dependence in 
this model. Considering the fugacity of the polymers 
as a temperature-like variable, bulk phase coexistence 
between a phase rich in colloids (analogous to a liquid 
in the liquid-vapor transition) with colloid packing frac- 
tion J]]., and a phase poor in colloids (analogous to the 
vapor) with colloid packing fraction 77^ < occurs on 
the line z^°°^{zp). The latter is determined by the equa- 
tion fic/ksT = /i™°''(zp)/fcBr in the plane of variables 
(zp, Zc). The variable A/i = /Xc — /x^"™ thus plays a role 
analogous to that of the magnetic field H in the Ising 
model, and ?7'(zp) — (■^p) corresponds to 2msp{T) for 
— > 0+, with msp(T) = — (9/sing/9-ff )t the sponta- 
neous magnetization or order-parameter. 

While in the Ising model there is a symmetry with re- 
spect to the sign of H, no corresponding symmetry with 
respect to the sign of Afj, occurs in colloid-polymer mix- 
tures, of course. In the Ising model, the two branches 
±msp(T) of the binodal, corresponding to positive and 
negative fields, are equivalent; this is the particle-hole 
symmetry of the lattice gas. Obviously, no such symme- 
try exists in asymmetric mixtures. Therefore, the above 
identification of variables holds only to a first approxima- 
tion, and for a more precise discussion of critical phenom- 
ena one needs to consider "field mixing " effects ka. The 
temperature-like variable t should then be taken parallel 
to the curve z^°^^{zp) at the point z"''(zp = Zp''*), and 
the field-like variable H becomes a suitable linear com- 
bination of the variables and Zp — z™'* (s^, E3 ■ 

Even more intricate becomes the identification of the 
variable corresponding to Hi, the local magnetic field 
coupling to the spins at the surface of the Ising magnet. 
The case Hi = physically means "neutral walls" , in 
which case the walls prefer neither the "spin-up" -phase 
nor the "spin-down" -phase of the ferromagnet. It is not 
at all obvious which choice of surface interactions would 
correspond to such a "neutral wall" , preferring neither 
the colloid-rich phase, nor the polymer-rich phase of our 
model. Of course, a completely analogous difficulty oc- 
curs in studies of capillary condensation and/or surface 
critical phenomena associated with the vapor-liquid criti- 
cal point of ordinary fluids 0,0,^^. We conclude that, 
in general, whatever the choice for the strength and range 
of the particle-wall interaction, it will likely correspond 



to a situation Hi ^ 0, but we have no possibility to pre- 
dict the strength (and even the sign) of Hi beforehand. 
In principle, a careful analysis of order parameter profiles 
at the bulk critical point of a very "thick" film (which ap- 
proximates a semi-infinite system with two surfaces) for 
various wall-particle interactions could provide insight 
into how to realize a situation with Hi = at criticality 
|38j . This, however, gives no guarantee that, for the same 
choice of interactions, one also has Hi = outside of the 
critical region. 

On the other hand, if our simulations display only very 
small shifts of fj,'^°°^{zp)/kBT with decreasing film thick- 
ness, one may assume that Hi in our model is indeed 
small. In this case, Eqs. ^ and lO still hold. However, 
if the data indicate that this is not the case, and rather 
the inverse limit HiD^^/'^ ^ 1 is reached instead, simple 
power laws would again result 

Uift = -D-^'''X,{^), AH, = -D-^/'Y^i^). (4) 

While the first power law of Eq. Q is the same as be- 
fore (only the constant of proportionality has changed) , 
the second power law clearly implies a somewhat faster 
decay than found in Eq. Q. Finally, for HiD^^/" k, 1 
a crossover behavior should be detectable, but to clearly 
identify such behavior may be difficult. In any case, the 
above discussion already shows that (confined) colloid- 
polymer mixtures may give rise to a much more com- 
plex behavior than the simple Ising model [Tll | , and hence 
their study should be rewarding. 



III. MODEL AND SIMULATION TECHNIQUE 

In this work, the colloids are modeled as hard spheres 
with diameter CTc, while the polymers are described as 
spheres with diameter Cp. The polymers may not over- 
lap with colloids, but there is no interaction between the 
polymers. The choice of this model motivated 
by the fact that flexible polymers form random coils with 
a rather large gyration radius i?g, which may interpene- 
trate at very low energy cost (in particular if the solution 
in which the polymers are dissolved, and the colloidal 
particles are suspended, is close to Theta-point condi- 
tions 01). Colloids, on the other hand, can be prepared 
with interactions of very short range [l^, ■ 

We specialize to a size ratio q = Upjuc = 0.8 in the fol- 
lowing and choose CTc = 1 as our unit of length. The poly- 
mer and colloid packing fractions arc rjp — napNp/ {QV) 
and ?7c = 'ko'^Nc/ (&V), when we have Np polymers and 
Nc colloids in our system of volume V , respectively. We 
choose a volume in the geometry of a rectangular box of 
linear dimensions x Ly x with = Ly ~ L and 
L2 = D, respectively. Periodic boundary conditions are 
applied in the x and y-directions. In the remaining z- 
dircction, we place two L x L walls, at z = and z = D, 
respectively, which act as hard walls for both polymers 
and colloids, while in addition a repulsive step potential 
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of height e (in units of fceT) acts on the colloids. The 
full colloid- wall interaction itcw(^) thus reads as: 

oo for h < 0, 
Mcw(/i) = "( e for < < crc/2, (5) 
otherwise, 

with h the distance from the surface of the colloid to 
the wall. In this work, e = 0, 0.25, 0.5, 1.0 and 2.0 are 
considered. Note that the case e = (purely hard walls 
for both colloids and polymers) and D = 5 has already 
been investigated in our previous work Q. In this case, a 
very strong attraction between the colloids and the walls 
develops due to the depletion effect 129 . 

Following common practice [T^ l29l l30l |. we choose 
the so-called (dimensionless) "polymer reservoir packing 
fraction" rj^ = TrzpCTp/G, rather than the polymer fugac- 
ity Zp, as the temperature- like variable. As in our pre- 
vious work \A I29I |30j , we apply a grand canonical clus- 
ter move |29j, together with a very efficient reweighting 
scheme, successive umbrella sampling |47| . in order to 
obtain the distribution function P^(77c|?7p, Zc), defined as 
the probability of observing the system with colloid pack- 
ing fraction rjc at "inverse temperature" r]^ and colloid 
fugacity Zc- By the subscript L, we remind the reader 
that quite generally there will be finite size effects, and a 
suitable extrapolation to the limit L — > cx), keeping the 
film thickness D fixed, may be required. While for states 
that are far away from phase coexistence, Pl(77c|'7pj ^c) 
is a single-peaked function, near two-phase coexistence 
a double-peak structure develops. The precise location 
of the fugacity Zc at which two-phase coexistence occurs 
is obtained using the equal-weight rule . The po- 

sitions of the two peaks in Pi(?7cb?p, 2;c) then yield esti- 
mates for the two branches, 77J. and 77^, of the binodal. In 
addition, we also study reduced moments of Pi(?7c|77p, Zc) 
at phase coexistence, defining an analogue of the order 
parameter of the Ising model. 



m = T]c- (/7c), iVc) = J ■q^PLiVcWp, Zc)dric (6) 
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FIG. 1: Moment ratio f/4 for a film of thickness D = 5 plotted 
versus the polymer reservoir packing fraction jf^ and several 
choices of the lateral dimension L = 12.5, 15.0, 17.5, 20.0, 
22.5, and 25.0. The intersection point yields an estimate for 
the critical "inverse temperature" »7p,cr = 0.892. The inset 
shows a log-log plot of the cumulant slope Y\ versus the linear 
dimension L. The straight line yields the exponent I/^'ch, 
from which we deduce v^^ — 0.933. 
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TABLE I: Critical parameters of the AO model for different 
film thicknesses D, compared to the corresponding 2d and 
3d bulk values. Note that the simulation data reported in 
this table were obtained using e — 0.5 as step-height in the 
colloid-wall interaction of Eq. 



coexistence near criticality is useful, and can yield even 
more accurate information on the critical behavior, the 
computer time resources for such a study in the present 
case are very demanding, and hence not attempted here. 



As is well known, following the behavior of moment ratios 
such as 



(8) 



along the path where Zc is at phase coexistence in the 
(zc, 77p)-plane, is useful for locating the critical point of 
the system; this is the so-called "cumulant intersection 
method" [53, ISll J52| . Although it has been demon- 
strated recently [a|53 that a study of the full distribution 
Pi(?7c|77p, Zc) and its moments also for fugacities Zc off 



IV. RESULTS FOR e = 0.5 

In this section, we present results using step-height 
e = 0.5 in the colloid-wall interaction of Eq. (0). We 
ultimately aim to test the Fisher-Nakanishi scaling pre- 
dictions |5( for the shift in the critical point parameters. 
To this end, the location of the critical point as function 
of the film thickness D, the binodal. and the behavior of 
the critical colloid packing fraction are investigated first. 
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FIG. 2: Same as Fig.0 but for film thickness D = 10, yielding 
the estimates 7?p,cr = 0.810 and i^off ~ 0.74. 




FIG. 3: Binodals of the AO-model with q = 0.8 in bulk (fuU 
curve without data points) and in confinement by parallel 
hard walls, to which a repulsive potential on the colloids with 
strength e = 0.5 was added (see Eq. JSJ). Circles denote data 
for D — 3, L — 18; squares for D — 5, L — 20; diamonds 
for L> = 75, L = 30; and triangles for D = 10, L = 30. The 
dotted lines denote the estimated coexistence diameters. 



A. Critical point 

We aim to locate the critical value of if^, for a number 
of thicknesses D. To this end, the cumulant ratio U^, 
sec Eq. ©, was measured as function of r^p using several 
lateral dimensions L. The results are shown in Figs. ^ 
and 12 for 13 = 5 and D = 10, respectively. We found 
that for D < 10, simulation data of meaningful accuracy 
could still be generated, although for D = 10 there is 
the need to choose the lateral linear dimension L as large 
as L = 40. From the cumulant intersection points, rj^ ^.j. 
is obtained. The "effective" correlation length exponent 
VcS is obtained from the L-dependence of the cumulant 
slope Yi evaluated at the intersection point. One expects 



that Yi cx 1/-^/'"="; the insets in Figs. and [3 show that 
our simulation estimates of Yi are indeed compatible with 
this relation. Additional simulations were performed for 
thicknesses D = 3 and D = 7.5; the resulting estimates 
of and i^cff , as well as the cumulant value t/| at the 
critical point, and the corresponding coexistence colloid 
chemical potential, arc collected in Table From these 
data, it is clear that, by increasing the film thickness D, 
a gradual crossover in the effective critical behavior from 
2d Ising toward 3d Ising universality occurs. While for 
D = 3, t'cff is very close to the exact value v = 1 of the 
2d Ising model, and similarly for the corresponding cu- 
mulant value ?7| (for the 2d Ising model, the very precise 
estimate ~ 0.856 was obtained H^l), with increasing 
D a clear shift toward 3d Ising values is observed. 

Of course, the smooth decrease of these "effective" 
values can only be taken as a very rough description 
of the theoretically expected crossover scaling scenario 
H i, m [13. For the case of D = 5 and hard- wall 
boundary conditions, a more elaborate analysis is pre- 
sented in Ref. . While in this case the cumulant slopes 
yield the same "effective" exponent i^eS ~ 0.933 as found 
here for D = 5 0, the more elaborate analysis in fact 
revealed Q that the asymptotic critical behavior is in 
reality described precisely by 2d critical exponents, but 
this can be seen only in a very narrow region around 
the critical point. Analyzing any such critical quantity 
on a log-log plot, one typically finds a systematic curva- 
ture: the slope of the curve approaches the 2d value for 
very small \t\ and then very gradually bends over in the 
direction towards the 3d value (the latter is not really 
observed because the crossover is not yet complete as |i| 
becomes so large that one leaves the critical region) p. 
Since a thorough analysis of crossover scaling requires an 
enormous investment of computer resources, we have re- 
frained from doing so, since there is no reason to expect 
any significant surprises. 



B. Binodal and coexistence diameter 

Next, we consider the binodal. Fig. 13 shows some of 
our "raw" data for the binodals in confinement, together 
with our estimates of the corresponding coexistence di- 
ameters (broken lines). By "raw" we mean that no finite 
size scaling analysis was performed on the data. The 
coexistence diameter S is defined as the average of the 
coexisting phase densities 



(9) 



where, as before, 77^ and ry' represent the coexistence col- 
loid packing fractions on the vapor and liquid branches 
of the binodal, respectively. The diameters terminate at 
the critical value of r]^, for which we used the estimates 
listed in Table HJ Since the data were obtained for finite 
L < 30, the two branches of the binodal do not merge at 
the critical point, but rather extend beyond it, bending 
over into the one-phase-region. This effect of "finite size 
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FIG. 4: Finite size coexistence diameters Sl,ct as function of 
for various choices of D as indicated. The arrow marks 
the diameter at the critical point in the thermodynamic limit 
of the bulk system, i.e. in the absence of walls, and was taken 
from Ref. 
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FIG. 5: Variation of the diameter at the critical point 
5* = lim_L^oo i^L.cr as function of D^^ . The inset shows the 
diameter 5 away from the critical point, choosing rf^ = 1.1, as 
function of . 



tails" or "rounding" of the order parameter into the dis- 
ordered region of the phase dia gram is well-known from 
simulations of the Ising model |50l Isil Is^ l . It is due to 
the fact that the order parameter distribution is clearly 
double peaked right at criticality and also in the one- 
phase region, if the lateral linear dimension L does not 
yet exceed the correlation length [s^, |^ 113 • The ex- 
ception in Fig. O is the full curve, which represents the 
binodal of the 3d bulk AO model. This curve was ob- 
tained using the finite size scaling approach of Ref. |53l| , 
and; on the scale of Fig.O should rather precisely reflect 
the true thermodynamic limit form L — > oo of the bulk 
binodal. 

In order to more accurately determine the colloid pack- 
ing fraction at the critical point 7?c,cn we have examined 
the i-dcpendcncc of 5l,ci, defined as the value of the co- 
existence diameter 5, given by Eq. @, as obtained in a 
finite system of lateral dimension L, at the critical value 
of r^p. The result is shown in Fig.^ where 6l,cv is plot- 
ted as function of 1/i, using for rf^ the values listed in 
Tabled The data display significant scatter, but rather 
intriguing behavior is revealed nevertheless. By increas- 
ing the film thickness _D, also 5l.ct: increases, away from 
the bulk value (arrow in Fig. ^J, although in the limit 
D oo precisely this bulk value should be recovered 
again. A possible explanation for this non-monotonic 
behavior with D, may be found in the precise critical be- 
havior of the coexistence diameter in the thermodynamic 
limit. The critical behavior of 6 in the limit L — > oo is 
rather intricate, and governed by several competing sin- 
gular terms [s^ 

S = 7,e,cr (1 + A2pt^'^ + A^^J^-" + A^t + . . ) , (10) 



colloid packing fraction at the critical point. Note that in 
2d, the first singular term has a much smaller exponent 
(2/3 = 1/4), but also the amplitude A2P may be very 
small compared to the amplitudes of the other terms. 
The next order term may involve a logarithmic correc- 
tion [5^ ; recall that a = in the 2d Ising model implies a 
logarithmic divergence of the specific heat. It is expected 
that in a finite system, the singular terms t^^ oc ^~2/3/i/ 
and i^^" oc ^-(1-")/'', cross-over to correction terms, 
with ^ ultimately replaced by L. Hence, we expect for 
the quantity ^^.cr in Fig. ^ the scaling form 



= ri,.AD)\^ + M[AD/L)L-^P/^ 
+ ii_„(i?/L)L-(i-")/'' 
+ ii(Z?/L)L-i/'' + •■•]. 



(11) 



with relative distance from the critical point 
Vp/Vp cr ~ li (nonunivcrsal) amplitudes Ai, and r/c. 



the 



Here, we have assumed that the amplitudes Ai are func- 
tions of the aspect ratio D/ L; such a dependence can be 
motivated by finite size scaling arguments 0, ^|. In 
view of the rather large error bars in S^ cr in Fig. 01 and 
the complicated structure of Eq. Hll|) . we feel that an 
accurate extrapolation to obtain 77c,cr(£') is not possible. 
Similarly, already in Ref. 0, it was pointed out that a 
reliable estimate of the critical behavior of the diameter 
in confinement from simulation data is not yet feasible. 
Hence, the only tentative conclusion we can draw from 
Fig. 2] is that presumably the dependence of rjc.cr{D) on 
D is non-monotonic: ijc,ci:{D) for small D does not differ 
much from 77c,cr(oo)j but with increasing D the difference 
increases first, reaches a maximum, and then decreases 
again. Such behavior could stem from competing signs 
of some of the amplitudes in Eqs. H1U|I and (|11|1 . 

For completeness, Fig. El shows 5* = limi^oo (^L.cr as 
function of -D~^, where S* was obtained by linearly ex- 
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FIG. 6; (a) Binodals of the AO-model with q = 0.8 in 
so-called grand-canonical representation, choosing the coexis- 
tence chemical potential of the colloids /ic""" and rfp as vari- 
ables. Shown is the bulk binodal (full curve), and the binodal 
in confinement for several film thicknesses D (broken curves) . 
The symbols mark the corresponding critical points, (b) Shift 
of the critical coexistence colloid chemical potential, given by 
Eq. (Uni, as function of ip-(^-^i)/''. (c) Shift of the criti- 
cal polymer reservoir packing fraction, given by Eq. as 
function of D'^^" . Linear straight line fits through the origin 
in (b) and (c) confirm that the data are compatible with the 
theoretical Fisher-Nakanishi predictions 1^. 



trapolating the data for Sl^ci- of Fig. 0] in the variable 
1/L. The resulting estimates of S* have also been col- 
lected in Tabled and reflect our best possible estimates 
of the critical colloid packing packing fraction 77c, cr in the 
thermodynamic limit (recall that Eq. H1U|I for the diam- 
eter reduces to r/c.cr = at the critical point t ^ 0). It 
is plausible from these data that the variation of S* with 
is indeed non-monotonic. However, the large error 
bars in d* , resulting from the uncertainty in the extrap- 
olation of Sl^ci to L — s- 00, prevent us from making any 
quantitative statements. In contrast, at values of rj^ that 
are far from the critical point for all choices of D con- 
sidered here, such as ij^ = 1.1, the L-dependence in the 
diameter can safely be neglected. In this case, the vari- 
ation with is clearly monotonic, and approximately 
linear for small D~^, see the inset of Fig. |S1 



C. Fisher-Nakanishi scaling and Kelvin equation 

Next, we consider the Fisher-Nakanishi scaling pre- 
dictions, essentially the main result of this paper. Re- 
call from the discussion in Section II, that the Fisher- 
Nakanishi scaling predictions pertain to the shift in the 
critical point parameters. More precisely, for the shift in 
the critical "inverse temperature" as function of the film 
thickness D, one expects that 

A^l,,{D) ^ f^l,,{D) rf^^,,{^) (X D-^/'^, (12) 



with u the correlation length critical exponent of the 
3d Ising model, rj^ ^^.{D) the critical value of if^ in the 
confined system of thickness D, and rj^ cr(c^) the criti- 
cal value of rjp in the bulk {D 00) system. The second 
scaling prediction of Fisher and Nakanishi pertains to the 
shift in the coexistence chemical potential of the colloids 
at criticality, as function of the film thickness D. It is 
expected that 

A^-f (i^) ^ /.™7(i?) - oc (13) 

with A and Ai the gap exponents introduced in Sec- 
tion II, i> again the correlation length exponent of the 
3d Ising model, fi^°^^{D) the coexistence chemical po- 
tential of the colloids at criticality in the confined sys- 
tem, and fil°°^ {00) the coexistence chemical potential of 
the colloids at criticality in the bulk system. Our re- 
sults have been collected in Fig. El The left panel shows 
the coexistence chemical potential ^1^°°^ of the colloids as 
function of r]^, for several values of the film thickness D. 
We remind the reader that the coexistence chemical po- 
tential of the colloids follows from the equal-weight-rule 
[isL I49I I . The shift in the coexistence chemical potential 
of the colloids at criticality, see Eq. (|13|l . is plotted in 
Fig. intb), as function of Here, the chemi- 

cal potential shifts were simply read-off from Fig. E^a); 
the corresponding values have also been collected in Ta- 
ble^ for completeness. The shift in the critical "inverse 
temperature", see Eq. (|12|l . is plotted in Fig. IH^c), as 
function of D"^^" , where for the critical values of rj^, 
the estimates of Table were used. In Figs. Elb) and 
(c), the known values of the exponents were used, as 
quoted previously. Validity of the Fisher-Nakanishi re- 
lations implies that the data collapse onto straight lines 
through the origin. Of course, data for small D, such 
as -D = 3 and D = 5, need not follow these relations, 
because the Fisher-Nakanishi predictions are valid only 
for asymptotically large D. For small Z?, corrections to 
finite size scaling likely come into play. In view of this, 
the agreement of our data with the scaling predictions, 
Eqs. (|2J) and |j2J), is already rather satisfactory. 

Finally, we consider the Kelvin equation, which is ex- 
pected to describe the shift in the coexistence colloid 
chemical potential away from the critical point. It is 
expected that 

(A^)Kelvin = l^r^'iD) - l^r^'ioo) OC 1/D, (14) 

with ^1°°^{D) the coexistence chemical potential of the 
colloids in the confined system, and /i™°''(oo) that in 
the bulk. Fig. [7| presents a test of the Kelvin equation, 
plotting the chemical potential difference as function of 
D~^, choosing three values of i]^ which are much higher 
than the corresponding critical values 77p (.^(£'). Once 
again, some systematic deviations at very small thick- 
nesses {D = 3 and Z) = 5) are seen, while for D = 7.5 and 
D = 10, one can already observe the expected straight 
lines through the origin. 
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FIG. 7: Test of the Kelvin equation. Shown is the chemi- 
cal potential difference, given by Eq. 1141 . as function of the 
inverse film thickness. Three values of 77p, chosen well above 
the critical values rjp^criD), are included. Broken straight lines 
show that the data are compatible with the Kelvin equation. 



D. Structural properties 

We now consider the local structure in the thin films. 
Figs. |S1 and |51 show density profiles of colloids and poly- 
mers (the latter are shown as insets), as a function of 
the distance z across the film. Results are shown for for 
D = 10 (Fig.© and D = 5 (Fig. III). The three values 
of rjp are chosen well away from the critical point, such 
that finite size effects in L can safely be neglected. The 
upper panel in each figure shows profiles obtained in the 
colloidal "vapor" phase; the lower panel in the colloidal 
"liquid" phase, that coexists with the vapor. It is clear 
that, despite the repulsive colloid-wall energy (e = 0.5 
in Eq. (Q), the entropically driven attraction of the col- 
loids to the walls still exists, and causes the formation of 
a rather dense layer of colloids at the walls. At not too 
large values of rj^, one can even see an indication of a 
second layer of colloids in the "vapor" phase. Of course, 
in the liquid phase, rather pronounced layering effects 
can always be detected. For the case D = 10, there are 
three colloidal layers adjacent to each wall; in the range 
3.5 < z < 6.5, the density profiles are essentially hori- 
zontal. This flatness indicates bulk-like behavior in the 
center of the thin fllm. In contrast, for the case D = 5, 
the layering associated with the left wall "interferes" with 
layering effects at the right wall, and no bulk region can 
be identified. For the particular parameters considered 
here. Fig. |51 shows that only five layers of colloids can be 
accommodated in a slit D = 5 colloid diameters wide. 
Note also the step-like increase in the colloid density at 
z = 1 and z = 4. This feature must be attributed to the 
step in the colloid- wall interaction, see Eq. jSJ. 

In view of the high colloid density reached locally in the 




FIG. 8; Colloid density profiles for film thickness D = 10 and 
lateral dimension L — 30. Profiles were obtained for three 
different values of r]p, as indicated, chosen to be well away 
from the critical point. Shown are profiles measured on the 
vapor branch of the binodal (a) and on the liquid branch (b) . 
The insets show the corresponding polymer density profiles. 



first layer adjacent to each wall, one may wonder whether 
this affects the lateral ordering of the colloids in planes 
parallel to the walls. Therefore, we have measured the ra- 
dial distribution functions gafsif) for colloid-colloid pairs 
(a/3 = cc) and colloid-polymer pairs (a/3 = cp), on the 
liquid branch of the binodal. Here, r denotes the magni- 
tude of the displacement vector r = \rij\ between pairs of 
particles = c,p). We have determined ga/sir) in the 
bulk region of the film (center), and in the vicinity of the 
walls, i.e. at a distance of 1 cr from the walls. Note that, 
in a system with walls (or virtual boundaries that define, 
e.g., the center of the confined system), the definition of 
the radial distibution function is different from that in 
the bulk. In an isotropic system with periodic boundary 
conditions, gapir) is normalized by a phase factor 47rr^. 
However, for the case that boundaries are present, this 
factor has to be modified such that it is then determined 
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FIG. 9: Same as Fig.|Hl but for D = 5 and L = 20. 

only by that part of the surface Anr"^ of a sphere around 
a particle that fits into the system. Thus, the normaliza- 
tion factor depends on the distance of the particles from 
the walls in z direction. The result for gap{r) is illus- 
trated in Fig. ^1 for thickness D ~ 10, lateral dimen- 
sion L ~ 30, and r/p = 1.1. The important conclusion of 
Fig.^lis that the radial distribution functions at the wall 
are almost indistinguishable from those of the bulk! This 
shows that, at the value of rf^ considered here, there is 
still no sign of wall-induced crystallization. Correspond- 
ing data for other film thicknesses, D ~ i and D = 5, 
have also been collected, but show the same trend as in 
Fig. EH and are therefore not shown. 



V. INFLUENCE OF THE COLLOID- WALL 
INTERACTION 

The results so far were obtained using e ~ 0.5 in the 
colloid- wall interaction of Eq. (jSJ. In this section, the pa- 
rameter e itself is varied, and the corresponding changes 




r 



FIG. 10: Radial distribution functions gcc{r) (main frame) 
and gcp{r) (inset) obtained in planes parallel to the confining 
walls (see details in text) for a film with thickness D — 10. 
The full curves show data obtained in the direct vicinity of 
the walls; the broken curves show data obtained in the bulk 
region of the film. Here, the bulk region was taken to be in 
the range 3.5 < z < 6.5, where the corresponding density 
profile is essentially flat, see Fig. |S1 



in the phase behavior will be discussed. Fig. ^2 shows 
the binodal for thickness D = 10 and D = 5, for two 
different values of e, compared to the bulk binodal. As 
before, e = implies pure hard walls We observe 
that, for large rj^, the binodal corresponding to e = 2.0 
and D = 10, approaches the bulk binodal rather closely. 
This indicates that e = 2.0 approximately cancels the en- 
tropic colloid-wall attraction. The latter is also evident 
from the corresponding density profiles, see Figs. I12r a) 
and (b). For e ~ 2.0, the peak in the colloid density pro- 
file in the layers adjacent to the walls, is significantly re- 
duced. In the liquid phase, we instead recognize a rather 
strong density jump for e ~ 2.0, at both z = 1 and 
z = D — 1, due to the step in the colloid- wall potential. 
Only for e = 0, do the density profiles remain smooth. 

Clearly, in future work, a careful analysis of more re- 
alistic colloid-wall interactions is desirable. Such a study 
is moreover warranted because we find, see Fig. El that 
variations in e significantly influence the radial distribu- 
tion functions, gccir) and 5cp(^), measured at the walls. 
By increasing e, a reduction of the first peak in gccir) oc- 
curs, with the second peak becoming almost completely 
washed out; for gcpif) the effect is even more dramatic. 
As expected, increasing e reduces the colloid density at 
the walls somewhat, thereby also weakening their ten- 
dency to order in the layers adjacent to the walls. Similar 
conclusions emerge from other wall thicknesses, too, and 
arc therefore not shown here. 
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FIG. 11: Binodal of the AO-model with q = 0.8 in bulk 
(full curve) and confinement (symbols). Shown are data for 
two system sizes; D = 10, L = 30 (a) and D = 5, L = 30 (b). 
For each system size, two values of the colloid-wall interaction 
parameter e, see Eq. are considered. Broken lines indicate 
the coexistence diameters of the confined systems. 



VI. DISCUSSION AND CONCLUDING 
REMARKS 

Colloid-polymer mixtures are important model sys- 
tems for the experimental exploration of fluid-fluid phase 
separation [l9ll20l | and interfacial phenomena 0,|22|. In 
order to guide possible future experiments on phase sep- 
aration of such systems in confined geometry, we have 
carried out a simulation study of a generic model sys- 
tem, the AO model, confined to a slit formed by two par- 
allel walls. The walls were not only impenetrable to both 
colloids and polymers, but in addition a short-ranged re- 
pulsive interaction between the walls and the colloids was 
added. Such an interaction can be realized, for instance, 
by a suitable coating of the walls with a moderately dense 
polymer brush (with the same chemical composition of 
the polymers dissolved in the solution). In this way, the 




FIG. 12: Colloid density profiles obtained in thin films at 
77p = 1.1, for two values of the film thickness D, and several 
values of the colloid-wall parameter e as indicated. Frames 
(a) and (b) show profiles obtained for D = 10, on the vapor 
and liquid branch of the binodal, respectively. Frames (c) and 
(d) show the corresponding profiles for thickness D = 3. The 
insets represent density profiles of the polymers. 
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FIG. 13: Radial distribution functions gcc{r) (main frame) 
and <7cp('') (inset), obtained in the boundary layer adjacent 
directly to the walls, for two values of e as indicated. Results 
were obtained for thickness D — 10, L = 30, and T^p = 1.1. 



strength of the depletion interaction between the colloid 
and the walls can be modified. 

Our results show that the unmixing binodal, describing 
phase separation in colloid-polymer mixtures, changes 
dramatically from its bulk form, when the width of the 
pore is in the range of 3-10 colloid diameters. While such 
ultra-thin pores seem to be experimentally realizable for 
colloidal systems, where the particle diameters are in the 
micrometer range, a corresponding study of a mixture 
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of small molecules or atoms, clearly, would be extremely 
difficult, if not impossible. Surprisingly, even for thin 
pores, we find that bulk-like phase behavior in the cen- 
ter of the pore is obtained, provided one is far away from 
criticality. Only for pores with thicknesses smaller than 
D = 3(Tc or so, can a bulk-like region in the center of the 
pore no longer be identified. In this case, the layering 
induced in the density by the left wall interferes with the 
layering due to the right wall. Therefore, such extremely 
thin slit pores cannot be accounted for quantitatively in 
terms of a decomposition of their properties into bulk and 
surface properties. Hence, the description of the shift of 
the transition between the phases poor in colloids ("va- 
por" ) and rich in colloids ( "liquid" ) in terms of the Kelvin 
equation, no longer is accurate. Nevertheless, agreement 
with the Kelvin equation is already recovered again above 
D ~ lOcTc or so, which is still remarkably thin. 

The main emphasis of this work, however, has been a 
study of the shift of the critical point. It has been shown 
that the scaling theory of Fisher and Nakanishi Q , orig- 
inally formulated for the Ising lattice gas model exhibit- 
ing full particle-hole symmetry between the coexisting 
phases, is compatible with our present numerical results. 
A nontrivial finding, which needs further investigation, is 
the fact that the critical value of the coexistence diame- 
ter approaches its limiting value in a non-monotonic way 



when D oo. 

In our study, we have obtained not only accurate bin- 
odal curves of confined colloid-polymer mixtures, but we 
have also characterized the coexisting colloid-rich and 
colloid-poor phases in terms of their concentration pro- 
files across the slit pore, and in terms of their radial 
distribution functions obtained in planes parallel to the 
walls. It would be interesting to extract corresponding 
information from experiments. Since the AO model is a 
crude simplification of reality, characteristic differences 
between simulation and experiment are to be expected, 
elucidating the limitations of the AO model. 

In future work, we plan to extend our study to con- 
finement between two parallel planar walls which are 
nonequivalent, whereby, for example, only one wall is 
coated by the polymer brush. In this way, an experimen- 
tal realization of the interface localization-delocalization 
transition Q may become possible. Also an extension of 
such studies to the dynamics of phase separation under 
confinement is of interest. We hope that the present work 
will stimulate corresponding experimental efforts. 
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